!==========Programa que calcula um pendulo forçado amortecido===========

PROGRAM exerB1

!=======================================================================

IMPLICIT NONE

!=======================Declaração de variáveis=========================

      INTEGER :: I, steps, flag
      REAL(8) :: tethai, tethai1, omegai, omegai1, Dt, currt
      REAL(8) :: l, totalt, Omega, gam, F0, periodo
      REAL(8), PARAMETER :: g = 9.80665d0, PI = 3.14159265358979d0

!=============Inicialização, leitura e abertura de arquivos=============

      WRITE(*,*) "Entre com o comprimento da haste:"
      READ(*,*) l
      WRITE(*,*) "Entre com o intervalo entre iterações delta_t:"
      READ(*,*) Dt
      WRITE(*,*) "Entre com o ângulo inicial do pêndulo (em graus):"
      READ(*,*) tethai
      WRITE(*,*) "Entre com o tempo total de simulação:"
      READ(*,*) totalt
      WRITE(*,*) "Entre com o parametro resistivo gamma:"
      READ(*,*) gam
      WRITE(*,*) "Entre com a amplitude da força externa F_0:"
      READ(*,*) F0
      WRITE(*,*) "Entre com a frequência da força externa:"
      READ(*,*) Omega

      !Abre arquivo para tetha(t)
      OPEN(UNIT=10, FILE='exerB1_out.dat', STATUS='unknown')

      !Abre arquivo para omega(t)
      OPEN(UNIT=20, FILE='exerB1_fout.dat', STATUS='unknown')

      !Converte o angulo de graus para radianos
      tethai = tethai*PI/180.0d0

      !Calculo do período
      periodo = (2.0d0*PI*DSQRT(l/g))*(1.0d0+(tethai**2.0d0)/16.0d0)

      tethai1 = 0.0d0
      omegai = 0.0d0
      omegai1 = 0.0d0
      currt = 0.0d0
      flag = 0

!=============================Faz as iterações==========================

      steps = totalt/Dt
      DO I=1,steps
            !Escreve no arquivo
            WRITE(10,*) currt, tethai
            WRITE(20,*) currt, omegai
            !Faz a iteração de omega
            omegai1 = omegai - (g/l)*DSIN(tethai)*Dt - gam*omegai*Dt +&
            &F0*DSIN(Omega*currt)*Dt
            !Faz a iteração de tetha
            tethai1 = tethai + omegai1*Dt
            !Muda os ângulos para deixar entre -pi e pi
            IF (tethai1 .GE. PI) THEN
                  tethai1 = tethai1 - 2.0d0*PI
            ELSE IF (tethai1 .LE. -1.0d0*PI) THEN
                  tethai1 = tethai1 + 2.0d0*PI
            END IF
            !Calcula a energia
            !Calcula em qual tempo ocorreu essa iteração
            currt = currt + Dt
            !Prepara variavaies para a proxima iteração
            omegai = omegai1
            tethai = tethai1
      END DO

      WRITE(*,*) "O período T em segundos é:", periodo

!=======================================================================

END PROGRAM exerB1

!=======================================================================
